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Abstract 


Nanoporous metals have been shown to exhibit radiation-tolerance due to the trapping of the 
defects by the surface. However, the behavior of vacancy clusters near the surface is not clear which 
involves the competition between the self-trapping and segregation of small vacancy clusters (V,) 
nearby the surface. In this study, we investigated the energetic and kinetic properties of small vacancy 
clusters near tungsten (0 0 1) surface by combining molecular statics (MS) calculations and object 
Kinetic Monte Carlo (OKMC) simulations. Results show that vacancies could be clustered with the 
reduced formation energy and migration energy of the single vacancy around a cluster as the 
respective energetic and kinetic driving forces. The small cluster has a migration energy barrier 


comparable to that for the single vacancy; the migration energy barriers for V -s and V7 are 1.80, 1.94, 


V. V. 


2.17, 2.78, 3.12 and 3.11 eV, respectively. Clusters ?and 3 become unstable near surface (0 0 1) 
and tend to dissociate into the surface. At the operation temperature of 1000 K, the single vacancy, V2, 
V; and V, were observed to segregate to the surface within a time of one hour. Meanwhile, larger 
clusters survived near the surface, which could serve as nucleating center for voids near the surface. 
Our results suggest that under a low radiation dose, surface (0 0 1) could act as a sink for small 
vacancy clusters, alleviating defect accumulation in the material under a low radiation dose. We also 
obtained several empirical expressions for the vacancy cluster formation energy, binding energy, and 


influence radius as a function of the number of vacancies in the cluster. 
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1. Introduction 

In future fusion nuclear reactors and advanced fission reactors, nuclear materials are exposed to 
a hostile environment involving intense heat fluxes, high temperatures and extreme fluxes of high- 
energy neutrons radiation [1]. Radiation exposure displaces atoms violently from their lattice sites, 
and produces point defects (interstitials and vacancies) and defect clusters, which may subsequently 
evolve into dislocation loops, voids and bubbles. Such microstructural change not only deteriorates 
the macroscopic properties of irradiated materials, but also severely degrades their performance and 
lifetime limits. In particular, voids, as a result of vacancy accumulation and growth, lead to swelling 
and eventual material failure [1-5]. Irradiation resistance is thus a major requirement for designing 
advanced nuclear materials in fusion reactors [6]. 

Concerning radiation tolerance, it is widely accepted that a variety of interfaces, e.g. surfaces in 
nanoporous (NP) materials and grain boundaries in nano-crystals can provide sinks for defects from 
radiation, and thus enhance the radiation tolerance [7-18]. For instance, irradiation experiments 
suggest that NP metals exhibit radiation-resistance [7-10]. By using in situ Kr-ion irradiation 
technique, Sun et al. observed frequent removal of various types of defect clusters by free surfaces in 
NP silver (Ag) [13]. They found, opposite to conventional wisdom, both types of diffusivities are 
lower in NP Ag due to the reduced interaction energy between isolated defect clusters in NP Ag. The 
void formation is expected to be suppressed as vacancies are absorbed by surfaces, since voids are 
developed through vacancy agglomeration [3]. However, the theoretical investigation of the vacancy 
cluster behavior in NP metals is limited (most of efforts have been made to investigate the vacancy or 
vacancy cluster behavior in the bulk material [19—25]). The evolution of the vacancy near the surface 
involves the competition of the clustering and segregation of vacancy clusters near the surface. 

As illustrated in Fig. 1, on one hand, vacancies could diffuse towards the surface and get trapped 
there. On the other hand, vacancies could also be clustered on the way to the surface. Determination 
of kinetics of defect clusters near the surface is essential to understand and model their migration and 
clustering in irradiated materials. In this work, NP tungsten (W) was chosen as a model system 
considering its engineering application as plasma facing materials and divertor armor in fusion 
reactors [26-28]. We investigated the energetic and kinetic properties of vacancy cluster and its 
nearby mono-vacancies in W on the atomic scale by using molecular statics (MS) simulations, 
including formation energy, migration barrier, binding energy and influence radius. With the energetic 
and kinetic parameters as input, the diffusion of the vacancy cluster near W (0 0 1) surface at 1000 K 
at long timescale was then investigated by object Kinetic Monte Carlo method (OKMC) [29]. We 
show that, although vacancies could be clustered near the surface, small clusters still segregate to the 
surface due to high mobility comparable to that for the single vacancy. 

2: Computational method 


2.1. Interatomic potential and surface models for MS simulations 


The bond-order potential (BOP) [30] was employed to describe atomic interaction among W 
atoms. The potential well describes the energetics and kinetics of point defects in W bulk. In bulk W, 
formation energy is 3.75 eV for the single vacancy and 10.19 eV for the self-interstitial atom of <1 1 
1> crowdion. The migration barrier of the vacancy is 1.80 eV and 0.002 eV for the self-interstitial 
atom. The potential has been used for studying radiation effects in W [17]. 

The BOP includes three-body contributions and angular dependence and thus it is time- 
consuming to perform simulations by using this potential. To increase the computational efficiency, 
we built two types of models. A bulk W model containing 5,488 atoms was constructed with a size of 
14x14x14(or 44.39x44.39x44.39 A), where ay is the lattice constant of 3.171 A. Periodic boundary 
conditions were imposed in all directions. The model has been checked to be large enough to 
investigate the energetic and kinetic behavior of small vacancy clusters. 

A (0 0 1) surface model containing 2,000 atoms was constructed with a size of 10a, in each 
direction. Periodic boundary conditions were imposed in [1 0 0] and [0 1 0] direction. A vacuum layer 
with a thickness of 15A was added in the [0 0 1] direction. The relaxed surface model is presented in 
Fig. 1. The size in the direction normal to the surface is large enough to avoid the interaction of the 
vacancy cluster with the other surface on the opposite side of the surface. 

2.2. Definition and calculations of defect properties 


The formation energy of a vacancy cluster was defined as: 


» (1) 


whereandare the total energies of the relaxed model with and without a vacancy cluster, 
respectively;is the atomic cohesive energy (—8.89 eV) of a bulk lattice W; n is the number of 
vacancies in the cluster. The system was relaxed by employing the steepest descent method. 


The binding energy of a vacancy with a vacancy cluster was defined as: 


» (2) 


whereactually is the bulk vacancy formation energy (3.75 eV). A position value of the binding energy 
indicates that it is energetically favorable for a mono-vacancy to locate near the cluster. 

The migration of small vacancy clusters in the bulk and near the surface was investigated by the 
nudged elastic band (NEB) method [31], which determined the saddle point and the migration barrier. 
There are 30 image states on the initial migration path. 

2.3. Details for OKMC simulations of vacancy evolution near the surface 

The kinetics of vacancy clusters near the surface was simulated by a standard OKMC method 
[29]. This technique allows us to follow the evolution of all the vacancy-type defects considered here, 
to take into account their spatial correlation and to cover the appropriate time and length scales. The 
input data for the OKMC include the energetic and kinetic properties of the vacancy clusters and the 


interaction range, which were generated by MS calculations. In the OKMC model, defects were 


considered as objects characterized by their space coordinates, nature (V and Vn in this study), shape 
(spherical here), mobility and dissociation rate. The diffusion processes for mobile defects were 
treated in a specific and efficient way: the details of the atomic transport were ignored and the jumps 
of the mobile species were bunched into single diffusion events. 

We used cubic simulation cells with a side of 20 nm. The upper and lower surfaces of the cell are 
defect sinks. As we investigated the behavior of a single vacancy cluster, the cluster was put at the 
center of the cell. The system temperature was chosen as 1000 K, which is the typical operation 
temperature of W in a fusion reactor [27]. As we explored the clustering and segregation of the 
vacancies near the surface, multiple vacancies were randomly generated. The initial concentrations of 
single vacancies were given two values of 100 and 300 appm, respectively. 


3. Results and discussion 


3.1. MS investigation of the formation of small vacancy clusters in W bulk 

We begin with finding the stable configuration for a vacancy cluster. Strictly speaking, for a 
given size of the cluster, it should reconsider the energy of all the different configurations to find the 
most stable one. However, such procedure can be very time consuming computationally due to 
multiple configurations of a vacancy cluster. The stable configuration for a vacancy cluster V,, was 
obtained by minimizing the system energy with a vacancy introduced at a particular site near the V n1. 
Moreover, we only searched the configurations for small vacancy clusters (Va, n=1-7). The 
corresponding configurations for the vacancy clusters can be seen in Fig. 2(a). We found that 
most of the vacancy clusters configurations are consistent with that in Ref. [30]. Note that, 
previously [30], Ahlgren et al. proposed that if any of the eight nearest atoms are missing for a 
particular W atom, the calculated total potential energy increases. The minimum energy configuration 
for a vacancy cluster should thus be a configuration where the total number of nearest neighbors in the 
lattice is at its maximum. Total energy calculations are thus needed only for those cluster 
configurations that result in a small number of broken W bonds. Nevertheless, the lowest energy 
configuration for a cluster of five vacancies is different from the square pyramid shape in 
Ref. [30], but consistent with the recent density functional theory result [32]. The formation 
energy for the present configuration is 1.28 eV lower than that in Ref. [30]. Therefore, such 
criteria may only apply for certain vacancy clusters. 

Fig. 2(b) shows the formation energies of vacancy clusters in the bulk as a function of the cluster 
size n. For a large void,is proportional to the surface energy of the void. Assuming the void to be 
spherical, it can be deduced thathas the form. As suggested in the continuum models describing void 
formation and growth during crystal growth [33], here parameter A is an effective surface energy for a 
cluster containing n vacancies. For the present small vacancy clusters, they should deviate from the 


classical n”? energy scaling behavior. Fitting the formation energies for vacancy clusters in bulk W 


with the same power-law results in the following equation: 
= 4.1562n°""??, (3) 


According to Eq. 2, the binding energy of a vacancy with a cluster is approximated by: 


EY = 4.1562[n°"” Ana) 83,79. (4) 

To check the prediction ability of this equation, we compared the MS results of the binding 
energy obtained according to Eq. 2 with the predictions by Eq. 4, as shown in Fig. 3. Note that, for a 
small V, (n<10), the binding energy was derived from V, formation energy (Eq. 2), while for a large 
Va, e.g. Vso, the bind energy was obtained by calculating the vacancy formation energy around the 
cluster (the reduction in the formation energy was binding energy). It can be seen from Fig. 3, Eq. 4 
well describes the vacancy binding energy with the V, as a function of n, although Eq. 4 overestimates 
the binding energy of the V with V2. 

In addition to the binding energy, we also calculated the influence range of the V,. Within the 
influence region, the V formation energy derivates from a bulk by at least 0.1 eV. By fitting the range 
to n, we found the range R approximately depends on n linearly as. The influence range for V, and V; 
was overestimated by such formula, but the formula well describes the vacancy cluster size 
dependence of R for other cases. 

3.2. MS investigation of the diffusion of small vacancy clusters in W bulk 

With the stable configurations of the clusters in hand, next we investigated their migration 
mechanisms to obtain the corresponding energy barrier. The vacancy clusters were assumed to 
migrate by successive nearest-neighbor mono-vacancy jumps. This motion may require 
passing through metastable intermediate states (see for example, V2, V4, Vs and V7). The 
migration energies are size-dependent: 1.80, 1.94, 2.17, 2.78, 3.12 and 3.11 eV for n = 1, 2, 3, 
4 and 7 respectively (Fig. 4). For V« (Fig. 2(a)), we found there are no one or two nearest- 
neighbor mono-vacancy jumps for Ve to migrate from one stable configuration to another one 
due to a high-symmetric octahedral configuration. We emphasize that the low values are obtained 
for V; and V4 comparable to that for the single vacancy. Larger V» clusters are believed to be 


immobile. 


To well understand the mobility of the small V, at an elevated temperature, we calculated the 
transition time of V, at 1000 K, which is a typical operation temperature for W in a fusion reactors 
[27]. As shown in Fig. 5, the time for small (n<5) to jump one step is less than 10’s, implying such V, 
has a not so low mobility. In other words,(n<5) can play an important role in cluster aggregation and 
growth or can diffuse into the surface. The transition time for V; and V; is longer than 10° s. Thus, the 
two clusters basically remain immobile at 1000 K. 

3.3. MS investigation of the formation and diffusion of mono-vacancies around the vacancy cluster 


To understand the clustering of the vacancy, we explored the interaction between a vacancy 


cluster and nearby mono-vacancies in terms of the influence of the vacancy cluster on the formation 
energy and migration energy barrier of nearby mono-vacancies. Fig. 6 presents the formation energies 
of mono-vacancies in the vicinity of the V, as a function of the distance of the vacancy from the 
cluster center. It can be seen, as a vacancy approaches to the cluster, its formation energy basically 
reduces. Thus, it is energetically favorable for a mono-vacancy to locate near a vacancy cluster and 
bind with the cluster into a new lager cluster. Meanwhile, the maximum reduction in the formation 
energy of mono-vacancies is cluster size dependent. 

The kinetics for mono-vacancies migration around vacancy clusters was further calculated by 
using the NEB method [31]. As shown in Fig. 7, the diffusion energy barrier for the single vacancy 
around near a V, depends on the migration path. There are always exist low energy barrier migration 
path for the vacancy to approach the V„. For example, the mono-vacancy has a barrier as low as 0.22 
eV on one path near V;. Such low energy barrier indicates the mono-vacancy far away from the 
cluster prefers to migrate towards the cluster kinetically. We also found that the energy barrier on 
some migration paths is higher than that for the single V diffusion in a perfect bulk W. Nevertheless, 
the average diffusion energy barrier for the V nearby the V, is reduced (Fig. 7(g)). Furthermore, it 
shows a decreasing trend with n, indicating that a larger V„ accelerates diffusion of the V more 
greatly. Meanwhile, the barrier for the V is lower than that for the V,. Thus, a Va could be considered 
as a Static sink for the V near itself. 

3.4. MS investigation of the diffusion of small vacancy clusters near W (0 0 1) surface 

A small V, (n=1-3) has a high mobility comparable to that for the single vacancy (Fig. 4). Next, 
we explore the diffusion and segregation of V and V; near a W (0 0 1) surface by calculating the 
formation energy and diffusion energy barrier of V and V3. Fig. 8 represents the results for V2. 
Initially, was placed at different sites near the surface as shown in Fig. 8(a). V2 contains two vacancies 
that are first-nearest neighbors. We designed several paths for V migration towards the surface. As V> 
is far away from the surface, its formation energy remains a constant bulk value. Only when V; is 
within a range of about 6.0 A, V, formation energy varies. Along the path AB-BC-CD-DE-EF-FG- 
GH-HI, V- first remains its bulk configuration and then becomes another configuration on the surface. 
The energy profile on this path is firstly increasing and then decreasing. The Vz formation energy with 
a configuration of FG is about 1.5 eV higher than that for a bulk V2. Along the path AB-BC-CD-DE- 
EF1-EH1-F1H1-G1HI1-HI1I1, V dissociates from a bulk configuration and then locates on the 
surface; vacancy clusters were also observed to dissociate near an iron grain boundary [34]. The 
energy profile along the path is basically decreasing. Therefore, it is therefore energetically favorable 
for V2 to locate at the surface. During the segregation process, Vz undergoes configurational 
transition. Furthermore, V2 on the surface contains two vacancies that are second-nearest neighbors, 
while in the bulk the two vacancies in V, are first-nearest neighbors. After a bulk Vz segregates to the 


atomic layers near the surface, the system energy decreases by about 3 eV. As V> transits from a 


configuration of GH (or G1H1) near the surface to a configuration of HI (or H1G1) on the surface, the 
system energy further reduces by about 3 eV. 

To show whether it is kinetically favorable for V» to segregate to the surface, we also calculated 
the migration energy landscape of V> near the surface. As shown in Fig. 8(c), V2 has a migration 
energy barrier of 1.94 eV as located in the bulk region (Fig. 4). But the energy landscape near the 
surface is rough. Except for the energy barrier from DF to EF/EF1, the barriers along other paths are 
smaller compared with the bulk diffusion energy barrier of V2. For example, V2 of FIH1 overcomes 
an energy barrier of 1.38 eV to locate at G1H1, which then nearly barrier-freely segregates to the 
surface along path G1H1-HI1I1. It is thus kinetically favorable for Vz to migrate towards the surface 
and reside therein. Previously [35], for the single vacancy, we have shown that the appearance of a 
high energy barrier near the surface causes the vacancy to try many times to segregate to the surface. 
We expect that aV, at AB will behave similarly, while a V2 at EF will take a directed walk towards the 
surface. 

Fig. 9 presents the energetic and kinetic properties of V; near the (0 0 1) surface. The initial 
configuration of V; is ABC, taken from Fig. 2. The computational procedure is similar to that for V2. 
We found both the formation energy profile and diffusion energy barrier landscape near the surface 
are similar to that for V2. After V; segregates from ABC to HHII on the surface along a diffusion- 
dissociation-segregation path, the system energy reduces by about 8.5 eV. Therefore, there is also an 
energetic driving force for V3 to migrate to the surface from the bulk. Near the surface, the energy 
barrier is significantly lower than that for V3 diffusion in the bulk. It is thus energetically and 
kinetically favorable for V; to enter the surface from the bulk. 

3.5. OKMC investigation of clustering and segregation of the vacancy near W (0 0 1) surface 

The above MS calculations provide atomic insight into the vacancy cluster behavior near the 
surface. To well understand the clustering and segregation of the vacancy near the surface at long 
timescale, we performed OKMC simulations at 1000 K. We found that, a single Vi, V2, V3 and V4, 
initially located at the center of a cubic cell with a side of 20 nm (in this case, the initial V, 
concentration is about 2 appm), all segregated to the surface within one hour, while Vs and V- took 
up to 10°s (about 10 days) to segregate to the surface. These results suggest that under a low radiation 
dose small vacancy clusters could not survive near the surface and will be absorbed by the surface at 
1000 K. 

To model the clustering and segregation of the V, under different levels of radiation dose, we 
designed another two OKMC simulations where the initial V concentration was assigned to be 100 
and 300 appm, respectively. After the system lasts for about one day, we examined the popular 
distribution of the vacancy clusters. We found that, as the initial concentration of the vacancy is 100 
appm, only several V; and Vs survived near the surface after one day. In the case of 300 appm, large 
vacancy clusters e.g. Vi. and Vi; appeared near the surface, implying that the clustering of vacancies 


and further growing of the cluster near the surface accompanied the segregation of small vacancy 
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clusters. Although the side of the cell is only 20 nm, vacancy clusters still accumulated near the 
surface. Therefore, surfaces in NP W may only efficiently alleviate defects accumulation in the 
system under a certain low radiation dose. This may be related to the small trapping range of the 
surface for the vacancy clusters (Figs. 8 and 9). 


4. Conclusion 


MS calculations and OKMC simulations have been combined to investigate the competition of 
the clustering and segregation of vacancies near the surface. Small vacancy clusters were found to 
have a high mobility comparable to that for the single vacancy. Mono-vacancies around the cluster 
prefer to bind with the cluster with reduced formation energy and diffusion energy barrier as diving 
forces. Clusterandbecome unstable near the (0 0 1) surface and get trapped by the surface via 
dissociation. At the operation temperature of 1000 K, Vi, V2, V3 and V, were observed to segregate to 
the surface within a time of one hour. Under a low radiation dose, surface (0 0 1) could act as a sink 
for small vacancy clusters, while under a high radiation dose, large vacancy clusters could survived 
near the surface and grow. We obtained several empirical expressions for the vacancy formation 
energy, binding energy and influence range as a function of the cluster size. Acknowledgments 
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Figure captions: 


Fig. 1. The (0 0 1) surface structure after relaxation. The position of the first atomic layer is indicated 
by the red dashed line at the top. The atoms are colored with their potential energies as indicated by 
the color bar on the right side. The vacancy (represented by a green cubic) near the surface could 


directly segregate to the surface or bind with vacancy clusters and then diffuse towards the surface. 


Fig. 2. (a) Stable configurations for small vacancy clusters. (b) Formation energies of the vacancy 


clusters in the bulk as a function of the number of vacancies n in the cluster. 


Fig. 3. Binding energy (left axis) and influence radius (right axis) for the vacancy cluster V, as a 


function of the number of vacancies, n. 


Fig. 4. Energy variation on the migration path for small vacancy clusters in W bulk. The 
configurations of the cluster before and after migration are given to show the migration mechanism. 
The vacant sites are represented by green squares. The normal lattice sites are marked by blue circles. 


The total energy of each system with a stable defect configuration is selected as the reference energy. 


Fig. 5. Migration barrier and transition time (left axis) for vacancy-type defects as a function of the 
number of vacancies n in the vacancy cluster. The transition time is calculated via where tọ with a 
value of 10°? s denotes the vibrational period of W atoms in the bulk. E, means migration energy 
barrier. The Boltzmann constant kz has a value of eV/K, and the absolute temperature T is set to 1000 


K, which is the typical operation temperature of W as a diverter. 


Fig. 6. Formation energy profile of the mono-vacancy near vacancy clusters as a function of the 
distance of the vacancy from the cluster center. The red dashed line indicates the bulk value of 
vacancy formation energy (3.75 eV). The black dashed lined indicates the minimal value of the 


formation energy. 


Fig. 7. (a) -(f) Migration of the mono-vacancy around a stable vacancy cluster (n=1, 2, 3, 5, 6, 7). The 
lowest energy configuration for the is indicated by green squares. A vacancy (indicated by a black 
cubic) nearby migrates towards the cluster and then binds with the cluster and becomes one 


component of the cluster (indicated by a red cubic). (g) The average diffusion energy barrier for the V 


near a V, as a function of n. 


Fig. 8. (a) The artificially designed path for the V migration towards the surface. The influence 


region of the surface is indicated by the two red dashed lines. (b) Formation energy profile for the 


cluster Vz near W (0 0 1) surface. (c) Diffusion energy landscapes for the V, migration towards the 


surface. The energy of the system with a bulk V- is chosen as reference. 


Fig. 9. (a) The artificially designed path for the V; migration towards the surface. The influence 
region of the surface is indicated by the two red dashed lines. (b) Formation energy profile for the 
cluster V; near W (0 0 1) surface. (c) Diffusion energy landscapes for the V; migration towards the 


surface. The energy of the system with a bulk V; is chosen as reference. 


Fig. 10. (a) Time needed for small vacancy cluster segregation to surface (0 0 1) at 1000 K drawn 
from OKMC simulations. (b) and (c) show the final population of vacancy clusters at 1000 K for an 
initial vacancy concentration of 100 and 300 appm, respectively. The evolution lasts for one day. In 
(b) and (c), vacancy clusters are shown by a green sphere with a radius proportional to the number of 


the vacancy in the cluster. 


